%% CLEAR ALL 
clear all;
close all;
clc;
rng(1, 'twister');
tic;

% ADD PATH 
% Determine where your m-file's folder is.
folder = fileparts(which("DGP_DFM.m")); 
% Add that folder plus all subfolders to the path.
addpath(genpath(folder)); 
% Saving path
file_dir = what('Results'); 
spath=file_dir.path;

% Set experiment
spec_id = 1;                % seed for random draws of specifications (= DGPs from encompassing model)
dgp_type = 'MP';            % structural shock: 'MP'
estimand_type = 'ObsShock'; % structural estimand: either 'ObsShock', or 'IV'
Interested_variable=21;     % 21=Industrial Production; 56=Unemployment rate
lag_type = 2;               % No. of lags to impose in estimation, or NaN (= AIC)
n_boot=500;                 % Number of bootstraps
weight=0.5;                 % Weighting for VAR IRFs
blocksize1=4;               % Size of block
blocksize2=20;              % Size of block
Ttest=100;                  % Size of test set (for cross-validation forecast)
% Settings
run('shared.m');
run(dgp_type);
run(estimand_type);


%% ENCOMPASSING DFM MODEL
% estimate DFM from dataset
DFM_estimate = DFM_est(DF_model.n_fac, DF_model.n_lags_fac, DF_model.n_lags_uar);

% extract and store estimated DFM parameters
DF_model.Phi           = DFM_estimate.Phi;
DF_model.Sigma_eta     = DFM_estimate.Sigma_eta;
DF_model.Lambda        = DFM_estimate.Lambda;
DF_model.delta         = DFM_estimate.delta;
DF_model.sigma_v       = DFM_estimate.sigma_v;
DF_model.variable_name = DFM_estimate.bplabvec_long;

% Represent as Model in ABCDEF Form
DF_model.n_s   = size(DF_model.Phi,2);
DF_model.n_eps = size(DF_model.Sigma_eta,2);
DF_model.n_y   = size(DF_model.Lambda,1);
DF_model.n_w   = size(DF_model.delta,1); % Warning: n_w stands for the dim of \omega_t
DF_model.n_e   = DF_model.n_w * DF_model.n_lags_uar;
DF_model.ABCD  = ABCD_fun_DFM(DF_model);

%% PREPARATION
% Select Individual DGPs from Encompassing Model
settings.specifications = pick_var_fn(DF_model, settings, spec_id);

% Compute True IRFs in Encompassing DFM
[DF_model.irf, settings.est.shock_weight] = compute_irfs(DF_model,settings);

% Compute Structural Target IRFs
switch estimand_type    
    case 'ObsShock'
        DF_model.normalized_irf = compute_normalized_irfs(DF_model,settings);
        DF_model.target_irf = DF_model.normalized_irf(settings.est.IRF_select, :);             
    case 'IV'
        DF_model.normalized_irf = compute_normalized_irfs(DF_model,settings);
        DF_model.target_irf = DF_model.normalized_irf(settings.est.IRF_select, :);
end

% Bayesian Local Projections Settings
nP=lag_type;
nH=settings.est.IRF_hor;
BLP_settings;

%Other settings
nlags=lag_type;
IRF_hor=settings.est.IRF_hor;
responseV = 3;
if strcmp(estimand_type, 'ObsShock')==1
normalizeV =1; %Obs
elseif strcmp(estimand_type, 'IV')==1
normalizeV =2; %IV
end






% Create matrices to store results
IRF_VAR_result=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_BLP_result=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_LP_result=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_result=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_fc_result=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);

IRF_VAR_UB_seqblock=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_VAR_LB_seqblock=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_LP_UB_seqblock=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_LP_LB_seqblock=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_UB_seqblock=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_LB_seqblock=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_UB_seqblock_bs1=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_LB_seqblock_bs1=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_UB_seqblock_bs2=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_LB_seqblock_bs2=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);

IRF_VAR_UB_bootrep=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_VAR_LB_bootrep=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_LP_UB_bootrep=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_LP_LB_bootrep=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_UB_bootrep=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_LB_bootrep=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);

IRF_VAR_seqblock_length=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_LP_seqblock_length=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_seqblock_length=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_seqblock_length_bs1=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_seqblock_length_bs2=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);

IRF_VAR_bootrep_length=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_LP_bootrep_length=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);
IRF_Pool_bootrep_length=zeros(settings.est.IRF_hor,settings.simul.n_MC,settings.specifications.random_n_spec);


correlation_varlp_seqblock_sim=zeros(2,2,settings.est.IRF_hor,settings.specifications.random_n_spec,settings.simul.n_MC);
correlation_varlp_seqblock_sim_bs1=zeros(2,2,settings.est.IRF_hor,settings.specifications.random_n_spec,settings.simul.n_MC);
correlation_varlp_seqblock_sim_bs2=zeros(2,2,settings.est.IRF_hor,settings.specifications.random_n_spec,settings.simul.n_MC);
correlation_varlp_bootrep_sim=zeros(2,2,settings.est.IRF_hor,settings.specifications.random_n_spec,settings.simul.n_MC);

%% Run Simulation 
parfor i_MC = 1:settings.simul.n_MC
    % Generate Data From Encompassing DFM
    rng(settings.simul.seed(i_MC), 'twister');
    data_sim_all = generate_data(DF_model,settings);

% Setting for Parfor (Local file)
IRF_VAR_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_BLP_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_LP_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_fc_spec=zeros(IRF_hor,settings.specifications.random_n_spec);

% Confidence bands Sequence Block Bootstrap
IRF_VAR_UB_seqblock_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_VAR_LB_seqblock_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_LP_UB_seqblock_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_LP_LB_seqblock_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_UB_seqblock_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_LB_seqblock_spec=zeros(IRF_hor,settings.specifications.random_n_spec);

% Confidence intervals Sequence Block Bootstrap (fix block size)
IRF_Pool_UB_seqblock_bs1_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_LB_seqblock_bs1_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_UB_seqblock_bs2_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_LB_seqblock_bs2_spec=zeros(IRF_hor,settings.specifications.random_n_spec);

% Confidence bands Bootstrap with replacement
IRF_VAR_UB_bootrep_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_VAR_LB_bootrep_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_LP_UB_bootrep_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_LP_LB_bootrep_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_UB_bootrep_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_LB_bootrep_spec=zeros(IRF_hor,settings.specifications.random_n_spec);

%Confidence interval length
IRF_VAR_seqblock_length_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_LP_seqblock_length_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_seqblock_length_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_seqblock_length_bs1_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_seqblock_length_bs2_spec=zeros(IRF_hor,settings.specifications.random_n_spec);

IRF_VAR_bootrep_length_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_LP_bootrep_length_spec=zeros(IRF_hor,settings.specifications.random_n_spec);
IRF_Pool_bootrep_length_spec=zeros(IRF_hor,settings.specifications.random_n_spec);

% Covariance
covariance_varlp_seqblock=zeros(2,2,IRF_hor);
covariance_varlp_seqblock_bs1=zeros(2,2,IRF_hor);
covariance_varlp_seqblock_bs2=zeros(2,2,IRF_hor);
covariance_varlp_bootrep=zeros(2,2,IRF_hor);

IRF_VAR_variance_seqblock=zeros(1,IRF_hor);
IRF_LP_variance_seqblock=zeros(1,IRF_hor);
IRF_VARLP_covariance_seqblock=zeros(1,IRF_hor);
IRF_VAR_variance_seqblock_bs1=zeros(1,IRF_hor);
IRF_LP_variance_seqblock_bs1=zeros(1,IRF_hor);
IRF_VARLP_covariance_seqblock_bs1=zeros(1,IRF_hor);
IRF_VAR_variance_seqblock_bs2=zeros(1,IRF_hor);
IRF_LP_variance_seqblock_bs2=zeros(1,IRF_hor);
IRF_VARLP_covariance_seqblock_bs2=zeros(1,IRF_hor);
IRF_VAR_variance_bootrep=zeros(1,IRF_hor);
IRF_LP_variance_bootrep=zeros(1,IRF_hor);
IRF_VARLP_covariance_bootrep=zeros(1,IRF_hor);

IRF_Pool_sd_seqblock=zeros(1,IRF_hor);
IRF_Pool_sd_seqblock_bs1=zeros(1,IRF_hor);
IRF_Pool_sd_seqblock_bs2=zeros(1,IRF_hor);
IRF_Pool_sd_bootrep=zeros(1,IRF_hor);

correlation_varlp_seqblock_spec=zeros(2,2,IRF_hor,settings.specifications.random_n_spec);
correlation_varlp_seqblock_spec_bs1=zeros(2,2,IRF_hor,settings.specifications.random_n_spec);
correlation_varlp_seqblock_spec_bs2=zeros(2,2,IRF_hor,settings.specifications.random_n_spec);
correlation_varlp_bootrep_spec=zeros(2,2,IRF_hor,settings.specifications.random_n_spec);

scale_VAR=zeros(1,1);
scale_LP=zeros(1,1);


% Setting end

    for i_spec = 1:settings.specifications.n_spec
       % Select data       
        data_sim_select = select_data_fn(data_sim_all,settings,i_spec);

% Common Setup for All Estimation Methods
% Prepare data
%run('Estimation_Setup'); % common setup for all estimation methods
if strcmp(estimand_type, 'ObsShock')==1
Y = [data_sim_select.data_shock,data_sim_select.data_y];
elseif strcmp(estimand_type, 'IV')==1
Y = [data_sim_select.data_z,data_sim_select.data_y];
end

N_Y=size(Y,2);
n_always=6;  
data_always=Y(:,1:n_always);

% Estimate IRFs by VAR
Y_var=data_always;
[Y_var_lagsnon1,Y_var_lags,Y_var_ad]=lag_variable(Y_var,nlags);
[T_var,N_var]=size(Y_var_ad);
[M_var]=size(Y_var_lags,2);

[IRF_VAR_full]=VAR_model(Y_var,IRF_hor-1,nlags,1,0);
IRF_VAR = IRF_VAR_full(responseV,:) / IRF_VAR_full(normalizeV,1); 
IRF_VAR_spec(:,i_spec)=IRF_VAR';

% Estimate IRFs by BLP
IRF_BLP_full=IRF_BLP_estimation(data_always,lag_type,IRF_hor,Y_var_lagsnon1,Y_var_lags,Y_var_ad,hyperPars,hyperPriorsOptions,0);
IRF_BLP = IRF_BLP_full(responseV,:) / IRF_BLP_full(normalizeV,1); 
IRF_BLP_spec(:,i_spec)=IRF_BLP';

% Estimate IRFs by LP
[IRF_LP_full]=LP_Jorda_model(Y_var,IRF_hor-1,nlags,1,0);
IRF_LP = IRF_LP_full(responseV,:) / IRF_LP_full(normalizeV,1); 
IRF_LP_spec(:,i_spec)=IRF_LP';

% Estimate IRFs by Pooling
IRF_Pool=weight*IRF_VAR+(1-weight)*IRF_LP;
IRF_Pool_spec(:,i_spec)=IRF_Pool';

% Weighting based on cross-validation forecast
[weightfcVAR, weightfcLP]=Weight_CVforecast(Y_var,nlags,Ttest,IRF_hor-1);

% Estimate IRFs by Pooling
IRF_Pool_fc=weightfcVAR(responseV,:).*IRF_VAR+weightfcLP(responseV,:).*IRF_LP;
IRF_Pool_fc_spec(:,i_spec)=IRF_Pool_fc';

% Sequence of block bootstrap
[IRF_VAR_boot_org,IRF_LP_boot_org,~,~]=VARLP_seqblockbootstrap(Y_var,nlags,IRF_hor,n_boot,responseV,normalizeV,0,0,0);

% Sequence of block bootstrap (fix-block size)
[IRF_VAR_boot_org_bs1,IRF_LP_boot_org_bs1,~,~]=VARLP_seqblockbootstrap(Y_var,nlags,IRF_hor,n_boot,responseV,normalizeV,0,1,blocksize1);
[IRF_VAR_boot_org_bs2,IRF_LP_boot_org_bs2,~,~]=VARLP_seqblockbootstrap(Y_var,nlags,IRF_hor,n_boot,responseV,normalizeV,0,1,blocksize2);

% Residual based bootstrap
[IRF_VAR_bootrep,IRF_LP_bootrep]=VAR_bootstrap_replacement(Y_var,IRF_hor-1,nlags,1,n_boot,responseV,normalizeV);

% Construct confidence intervals
scale_VAR=1/ IRF_VAR_full(normalizeV,1); 
scale_LP =1/ IRF_LP_full(normalizeV,1);
zvalue=1.96;


for j=1:IRF_hor
covariance_varlp_seqblock(:,:,j)=cov(IRF_VAR_boot_org(:,j),IRF_LP_boot_org(:,j));
covariance_varlp_seqblock_bs1(:,:,j)=cov(IRF_VAR_boot_org_bs1(:,j),IRF_LP_boot_org_bs1(:,j));
covariance_varlp_seqblock_bs2(:,:,j)=cov(IRF_VAR_boot_org_bs2(:,j),IRF_LP_boot_org_bs2(:,j));
covariance_varlp_bootrep(:,:,j)=cov(IRF_VAR_bootrep(:,j),IRF_LP_bootrep(:,j));

correlation_varlp_seqblock_spec(:,:,j,i_spec)=corrcoef(IRF_VAR_boot_org(:,j),IRF_LP_boot_org(:,j));
correlation_varlp_seqblock_spec_bs1(:,:,j,i_spec)=corrcoef(IRF_VAR_boot_org_bs1(:,j),IRF_LP_boot_org_bs1(:,j));
correlation_varlp_seqblock_spec_bs2(:,:,j,i_spec)=corrcoef(IRF_VAR_boot_org_bs2(:,j),IRF_LP_boot_org_bs2(:,j));
correlation_varlp_bootrep_spec(:,:,j,i_spec)=corrcoef(IRF_VAR_bootrep(:,j),IRF_LP_bootrep(:,j));

% For sequence of block bootstrap (flexible block size)
IRF_VAR_variance_seqblock(:,j)=covariance_varlp_seqblock(1,1,j);
IRF_LP_variance_seqblock(:,j)=covariance_varlp_seqblock(2,2,j);
IRF_VARLP_covariance_seqblock(:,j)=covariance_varlp_seqblock(2,1,j);

% For sequence of block bootstrap (fix-block size)
IRF_VAR_variance_seqblock_bs1(:,j)=covariance_varlp_seqblock_bs1(1,1,j);
IRF_LP_variance_seqblock_bs1(:,j)=covariance_varlp_seqblock_bs1(2,2,j);
IRF_VARLP_covariance_seqblock_bs1(:,j)=covariance_varlp_seqblock_bs1(2,1,j);

IRF_VAR_variance_seqblock_bs2(:,j)=covariance_varlp_seqblock_bs2(1,1,j);
IRF_LP_variance_seqblock_bs2(:,j)=covariance_varlp_seqblock_bs2(2,2,j);
IRF_VARLP_covariance_seqblock_bs2(:,j)=covariance_varlp_seqblock_bs2(2,1,j);

% For residual based bootstrap
IRF_VAR_variance_bootrep(:,j)=covariance_varlp_bootrep(1,1,j);
IRF_LP_variance_bootrep(:,j)=covariance_varlp_bootrep(2,2,j);
IRF_VARLP_covariance_bootrep(:,j)=covariance_varlp_bootrep(2,1,j);


%Variance of Pooled IRFs
IRF_Pool_sd_seqblock(:,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_seqblock(:,j)+(1-weight)^2*scale_LP^2*IRF_LP_variance_seqblock(:,j)...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_seqblock(:,j));

IRF_Pool_sd_seqblock_bs1(:,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_seqblock_bs1(:,j)+(1-weight)^2*scale_LP^2*IRF_LP_variance_seqblock_bs1(:,j)...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_seqblock_bs1(:,j));

IRF_Pool_sd_seqblock_bs2(:,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_seqblock_bs2(:,j)+(1-weight)^2*scale_LP^2*IRF_LP_variance_seqblock_bs2(:,j)...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_seqblock_bs2(:,j));

IRF_Pool_sd_bootrep(:,j)=sqrt(weight^2*scale_VAR^2*IRF_VAR_variance_bootrep(:,j)+(1-weight)^2*scale_LP^2*IRF_LP_variance_bootrep(:,j)...
                              +2*weight*(1-weight)*scale_VAR*scale_LP*IRF_VARLP_covariance_bootrep(:,j));

end

% Confidence bands Sequence Block Bootstrap
IRF_VAR_UB_seqblock_spec(:,i_spec)=(IRF_VAR+zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_seqblock))';
IRF_VAR_LB_seqblock_spec(:,i_spec)=(IRF_VAR-zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_seqblock))';

IRF_LP_UB_seqblock_spec(:,i_spec)=(IRF_LP+zvalue*sqrt(scale_LP^2*IRF_LP_variance_seqblock))';
IRF_LP_LB_seqblock_spec(:,i_spec)=(IRF_LP-zvalue*sqrt(scale_LP^2*IRF_LP_variance_seqblock))';

IRF_Pool_UB_seqblock_spec(:,i_spec)=(IRF_Pool+zvalue*IRF_Pool_sd_seqblock)';
IRF_Pool_LB_seqblock_spec(:,i_spec)=(IRF_Pool-zvalue*IRF_Pool_sd_seqblock)';


% Confidence intervals Sequence Block Bootstrap (fix block size)
IRF_Pool_UB_seqblock_bs1_spec(:,i_spec)=(IRF_Pool+zvalue*IRF_Pool_sd_seqblock_bs1)';
IRF_Pool_LB_seqblock_bs1_spec(:,i_spec)=(IRF_Pool-zvalue*IRF_Pool_sd_seqblock_bs1)';

IRF_Pool_UB_seqblock_bs2_spec(:,i_spec)=(IRF_Pool+zvalue*IRF_Pool_sd_seqblock_bs2)';
IRF_Pool_LB_seqblock_bs2_spec(:,i_spec)=(IRF_Pool-zvalue*IRF_Pool_sd_seqblock_bs2)';


% Confidence bands Bootstrap with replacement
IRF_VAR_UB_bootrep_spec(:,i_spec)=(IRF_VAR+zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_bootrep))';
IRF_VAR_LB_bootrep_spec(:,i_spec)=(IRF_VAR-zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_bootrep))';

IRF_LP_UB_bootrep_spec(:,i_spec)=(IRF_LP+zvalue*sqrt(scale_LP^2*IRF_LP_variance_bootrep))';
IRF_LP_LB_bootrep_spec(:,i_spec)=(IRF_LP-zvalue*sqrt(scale_LP^2*IRF_LP_variance_bootrep))';

IRF_Pool_UB_bootrep_spec(:,i_spec)=(IRF_Pool+zvalue*IRF_Pool_sd_bootrep)';
IRF_Pool_LB_bootrep_spec(:,i_spec)=(IRF_Pool-zvalue*IRF_Pool_sd_bootrep)';

%Confidence interval length
IRF_VAR_seqblock_length_spec(:,i_spec)=(2*zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_seqblock))';
IRF_LP_seqblock_length_spec(:,i_spec)=(2*zvalue*sqrt(scale_LP^2*IRF_LP_variance_seqblock))';
IRF_Pool_seqblock_length_spec(:,i_spec)=(2*zvalue*IRF_Pool_sd_seqblock)';
IRF_Pool_seqblock_length_bs1_spec(:,i_spec)=(2*zvalue*IRF_Pool_sd_seqblock_bs1)';
IRF_Pool_seqblock_length_bs2_spec(:,i_spec)=(2*zvalue*IRF_Pool_sd_seqblock_bs2)';


IRF_VAR_bootrep_length_spec(:,i_spec)=(2*zvalue*sqrt(scale_VAR^2*IRF_VAR_variance_bootrep))';
IRF_LP_bootrep_length_spec(:,i_spec)=(2*zvalue*sqrt(scale_LP^2*IRF_LP_variance_bootrep))';
IRF_Pool_bootrep_length_spec(:,i_spec)=(2*zvalue*IRF_Pool_sd_bootrep)';

    end


% Store results
% IRFs
IRF_VAR_result(:,i_MC,:)=IRF_VAR_spec;
IRF_BLP_result(:,i_MC,:)=IRF_BLP_spec;
IRF_LP_result(:,i_MC,:)=IRF_LP_spec;
IRF_Pool_result(:,i_MC,:)=IRF_Pool_spec;
IRF_Pool_fc_result(:,i_MC,:)=IRF_Pool_fc_spec;


% Confidence intervals
IRF_VAR_UB_seqblock(:,i_MC,:)=IRF_VAR_UB_seqblock_spec;
IRF_VAR_LB_seqblock(:,i_MC,:)=IRF_VAR_LB_seqblock_spec;
IRF_LP_UB_seqblock(:,i_MC,:)=IRF_LP_UB_seqblock_spec;
IRF_LP_LB_seqblock(:,i_MC,:)=IRF_LP_LB_seqblock_spec;
IRF_Pool_UB_seqblock(:,i_MC,:)=IRF_Pool_UB_seqblock_spec;
IRF_Pool_LB_seqblock(:,i_MC,:)=IRF_Pool_LB_seqblock_spec;

% Confidence intervals Sequence Block Bootstrap (fix block size)
IRF_Pool_UB_seqblock_bs1(:,i_MC,:)=IRF_Pool_UB_seqblock_bs1_spec;
IRF_Pool_LB_seqblock_bs1(:,i_MC,:)=IRF_Pool_LB_seqblock_bs1_spec;
IRF_Pool_UB_seqblock_bs2(:,i_MC,:)=IRF_Pool_UB_seqblock_bs2_spec;
IRF_Pool_LB_seqblock_bs2(:,i_MC,:)=IRF_Pool_LB_seqblock_bs2_spec;

% Confidence bands Bootstrap with replacement
IRF_VAR_UB_bootrep(:,i_MC,:)=IRF_VAR_UB_bootrep_spec;
IRF_VAR_LB_bootrep(:,i_MC,:)=IRF_VAR_LB_bootrep_spec;
IRF_LP_UB_bootrep(:,i_MC,:)=IRF_LP_UB_bootrep_spec;
IRF_LP_LB_bootrep(:,i_MC,:)=IRF_LP_LB_bootrep_spec;
IRF_Pool_UB_bootrep(:,i_MC,:)=IRF_Pool_UB_bootrep_spec;
IRF_Pool_LB_bootrep(:,i_MC,:)=IRF_Pool_LB_bootrep_spec;


%Confidence interval length
IRF_VAR_seqblock_length(:,i_MC,:)=IRF_VAR_seqblock_length_spec;
IRF_LP_seqblock_length(:,i_MC,:)=IRF_LP_seqblock_length_spec;
IRF_Pool_seqblock_length(:,i_MC,:)=IRF_Pool_seqblock_length_spec;
IRF_Pool_seqblock_length_bs1(:,i_MC,:)=IRF_Pool_seqblock_length_bs1_spec;
IRF_Pool_seqblock_length_bs2(:,i_MC,:)=IRF_Pool_seqblock_length_bs2_spec;

IRF_VAR_bootrep_length(:,i_MC,:)=IRF_VAR_bootrep_length_spec;
IRF_LP_bootrep_length(:,i_MC,:)=IRF_LP_bootrep_length_spec;
IRF_Pool_bootrep_length(:,i_MC,:)=IRF_Pool_bootrep_length_spec;


% Average correlation across specifications
correlation_varlp_seqblock_sim(:,:,:,:,i_MC)=correlation_varlp_seqblock_spec;
correlation_varlp_seqblock_sim_bs1(:,:,:,:,i_MC)=correlation_varlp_seqblock_spec_bs1;
correlation_varlp_seqblock_sim_bs2(:,:,:,:,i_MC)=correlation_varlp_seqblock_spec_bs2;
correlation_varlp_bootrep_sim(:,:,:,:,i_MC)=correlation_varlp_bootrep_spec;



end % End simulation


% Estimate all results
Run_DFMDGP_all_results

% Save results
Run_DFMDGP_save_results
